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Abstract. A new parametrization of the Eilenberger equations of superconductivity in 
terms of the solutions to a scalar differential equation of the Riccati type is introduced. 
It is shown that the quasiclassical propagator, and in particular the local density of 
states, may be reconstructed, without explicit knowledge of any eigenfunctions and 
eigenvalues, by solving a simple initial value problem for the linearized Bogoliubov-de 
Gennes- equations. The Riccati parametrisation of the quasiclassical propagator leads 
to a stable and fast numerical method to solve the Eilenberger equations. 

1 Introduction 

According to the BCS theory of superconductivity the quasiparticle excitations 
above the Cooper pairing groundstate depend on spin ( f or J, ) and also on a 
particle-hole index (+ or — ) which indicates the flight direction of a quasiparticle 
(parallel or antiparallel to the Fermi velocity ). Coherent superpositions of 
such excitations form wave packets that transport energy, momentum, charge 
and spin inside a superconductor. 

In metals and alloys of interest to technical applications of superconductiv- 
ity the Cooper pairs display an even parity symmetry (spin singlet), and often 
the influence of paramagnetic effects ( Zeeman splitting, Pauli limiting, spin- 
orbit coupling etc.) may be ignored. Then spin and particle-hole indices may be 
identified. As a result the 4 x 4-matrix equations of superconductivity may be 
simplified to 2 x 2-matrix equations. 

In the following we use notation such that r refers to a point in position space 
(center of mass of a Cooper pair) , and Pf = fikp denotes a point on the Fermi 
surface FS . 

It is known that the characteristic length to heal a local (static) perturbation 
of the Cooper pairing amplitude A(r, pp) in a superconductor (due to the pres- 
ence of an impurity, a vortex line, an interface etc. ) is approximately £ = ; 
and often the quasiclassical condition ^ 1 is fulfilled. 

Then, as first shown by Eilenberger[l] and Larkin and Ovchinnikov[2] , the 
relevant part of the physical information coded in quantum mechanical expecta- 
tion values ( for example the charge density, the current, the pressure functional 
etc.) may be calculated more efficiently with the help of the quasiclassical prop- 
agator 




(1) 
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The quasiclassical propagator is, by definition, just the Green's function of the 
Gorkov theory of superconductivity in a form where it has been integrated with 
respect to the kinetic energy of the quasiparticles[4]. Remarkably, g(r;pp,ie n ) 
may be also calculated directly solving a transport type system of ordinary 
differential equations (the right hand side is a commutator): 



ifiv F -Wg(r;p F ,is n ) = 
ie n + v F • f A (r) -A(r,p F ) 



g(r;p F ,ie n ) 



(2) 



z}t( r ,p F ) -ie n - v F • § A (r) 
The physical solution to this equation must also fulfill a normalisation condition: 
g(r; p F , ie n ) ■ g(r; p F ,ie n ) = -n 2 • 1 (3) 



General symmetries of the Gorkov Green's functions imply corresponding sym- 
metries of the quasiclassical propagator: 



f(r;p F ,is n ) = -f*(r;p F ,-ie n ) 
g(r; p F ,is n ) = g(r; -p F , -ie n ) 
f(r;-PF,-ie n ) = f{r;p F ,ie n ) 
g(r;p F ,is n ) = g*(r; p F , -ie n ) 



(4) 
(5) 
(6) 

(7) 



In equilibrium the quasiclassical propagator also displays a particle-hole symme- 
try: 

g(r; p F ,ie n ) = -g(r; p F , ie n ) (8) 

This means that the trace of g(r;p F ,ie n ) vanishes. As a traceless 2 x 2-matrix 
the square of g should be equal to a multiple of unity: 



g(r; p F ,ie n ) ■ g(r; p F ,ie n ) = C ■ 1 



(9) 



Using the fact, that g 2 is a solution to the Eilenberger equations (provided g 
is a solution) , it follows that — ihvp-VC = 0, i.e. the scalar C is necessarily a 
constant along a straight line orientated parallel to the Fermi velocity v F . But C 
could still be a function of the form C = C (rAv F ;p F , ie n ) . The normalisation 
condition Eq.(3) fixes C such that g 2 = —it 2 ■ 1 for all straight lines orientated 
parallel to v F , and this for all Fermi momenta p F on the Fermi surface and 
also for all Matsubara frequency ie n . The particular value C — — n 2 is chosen 
in order to achieve consistency with the functional form of the quasiclassical 
propagator in the bulk. 

In thermal equilibrium the pair potential A(r, p F ), the electrical current J(r) 
associated with a (stationary) flow of quasiparticles, the local density of states 
N(r,E), the Gibbs free energy Gs of the superconducting state for weak cou- 
pling^], and other observables may be directly calculated using the quasiclassical 
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propagator: 



A(r,p F )= [ dp' F N FS (p F )V(p Fl p F )-k B T V f(v,p' F ,is n ) (10) 

Jfs \e n \< Uc 

J(r) = ^E/ ^Pf^(Pf) v f s(r,p' F ,ie n ) (11) 
N(r,E) = -- [ dp' F N(p' F )Img(r,p' F ,ie n ^ E + i0+) (12) 

(13) 

-2T ■ J^dE N(r,E) -\n(e^ + 
G S {T) = J dv +J FS dp F J FS dp' F AHr,p F )o(V- 1 ) pF ^ o4(r,p' F ) 

+ll(VAA(r)-B elt (r)) 2 

In these expressions the function N F s(p F ) denotes the (angle resolved) density 
of states in the normal phase at the Fermi level. This function typically enters 
as a weight function into Fermi surface integrals ( FS denotes the Fermi surface) 
of the Eilenberger propagator. In the isotropic case N F s(p F ) simplfies to the 
usual constant N(0). 

The calculation of Fermi surface integrals of the Eilenberger propagator be- 
comes comparatively simple in the bulk, where the pair potential, A(p F ), is 
independent on position r, and where the quasiclassical propagator assumes the 
form: 



g(p F ,is n ) = 



— 7T 



^sl + \A( PF )\ 



_ ( ie n ~A(p F ) 

2 ^(Pf) 1 -i£n 



(14) 



A considerably more complicated problem is posed when the pair potential 
depends on position r, for instance near a surface, in the vicinity of an implanted 
impurity or ion, or around a flux line in a type-II superconductor. 

Usually the solution g(r;p F , ie n ) of the Eilenberger equations must be found 
numerically. But the task is more difficult then just solving a differential equa- 
tion. To determine the pair potential A(r, p F ) and the magnetic field B(r) = 
V A A(r) from the (magnetostatic) Maxwell Equation , V A B(r) = ^LJ(r), one 
needs to solve a (nonlinear) selfconsistency problem, since J(r) and A(r,p F ) 
depend themselves on g(r; p F , ie n ). 



2 Eilenberger Equations along a Characteristic Line 



First we consider a layered material (normal axis parallel to c ) assuming, for 
example, a Fermi velocity v F that is orientated predominantly within the ab— 
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plane (Fermi circle). Let the triade |a, b, c| span an ortho-normal basis in the 
lab frame, while 9 denotes the angle the Fermi velocity v F makes with the a-axis. 
Clearly, along a straight line 

r(a;) = xv+yu (15) 
= r a (x)a+r b (x)b (16) 

— OO < X < oo 

with v and u denoting unit vectors (orientated parallel and orthogonal to vp , 
respectively) , 

v = cos(0)a + sin(0)b (17) 
u = - sin(0)a + cos(0)b , (18) 

the directional derivative vp ■ V in the Eilcnberger Equation Eq.(2) is equivalent 
to an ordinary derivative: 

d 

hv F • V<7(r;p F ,ie„) = hv F —g[r(x);p F ,ie n ] (19) 

The ^-dependent parameter y associated with such a characteristic line r(x) 
(see Eq.(15)) has the natural meaning of an impact parameter. The straight 
line r(x) intersects with a fixed position point 

r = r a k + r b b (20) 

( there where the solution g(r;pF,ie n ) is sought) at the particular parameter 
value x = xp . Introducing polar coordinates, 

r a + in = \]rl+r 2 b e* (21) 

it is evident that 

r a (x) + ir b (x) = (x + iy)e t6 (22) 

and 



x P + iy = ^rl + r\ (23) 

The extension to 3-dimcnsions is straightforward. For instance, for a spherical 
Fermi surface the unit vectors v and u are parametrised by two angles, the 
azimutal angle 9 e [0, 27r) and the polar angle \ € [0,tt), respectively: 

v = sin (x) cos(6')a + sin(#)b + cos (x) c (24) 

-1 d _ 

u = sin (x) — sin(#)a + cos(0)b = — v (25) 
L J o9 

Again, along a straight line, r(x) — xv + yu+zv A u , the directional derivative 
in the Eilenberger equations becomes just an ordinary derivative. Making the 
identification r =r a & + r b h + r c c = r(xp) explicit expressions for xp and both 
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'impact' parameters y and z in terms of 9, \ and the cartesian coordinates r a , 
ri, , r c of the fixed point r in position space are easily derived. 

Finally we simplify our notation by dropping the functional dependence of 

A 

A , A and 9 on arguments that stay constant as x varies from — oo to oo : 

A(x) = A[r(x),p F ] (26) 
is n (x) = ie n + vf • -A [r(x)] (27) 
g(x) = g[r(x);p F ,ie n ] (28) 



3 Riccati Parametrisation of Eilenberger Propagator 



Any traceless 2 x 2-matrix may be expanded into the basis 

1_ 



K 3 
K± 



7 T 3 



-- • (fi ±i? 2 ) 



(ri, t 2 ,and % are standard 2 x 2-Pauli matrices). We note that 

K+,K-] =-2^ 3 

'k 3 ,k±\ =±k ± 



(29) 
(30) 

(31) 
(32) 



The Eilenberger equations may then be rewritten along a characteristic line r(x) 
orientated parallel to the Fermi velocity in the form: 



ox 



-2e n (x)K 3 + A(x)K+ - A\x)K_ , g(x) 



(33) 



Let us consider the following 2x2 system of ordinary differential equations 
for an auxiliary propagator Y(x) (fundamental system): 



hv F ^Y(x) - (-2s n {x)K 3 + A(x)K+ - A\x)K^j Y{x) 
?(0) - Y 



(34) 
(35) 



The initial values for Y(x) at i = may be prescribed in terms of a (yet 
unknown) constant 2x2 matrix Yq of rank 2 . We may reconstruct the physical 
propagator g , the one that solves the Eilenberger equations and respects the 
normalization condition, g(x) ■ g(x) = —it 2 ■ 1 , from the fundamental system 

Y(x): 

g{x) = -m ■ Y(x) ■ 2K 3 ■ Y'^x) (36) 

By putting x at the end of the calculations to the particular value xp , the phys- 
ical propagator (i.e. the input into the selfconsistency equations) is recovered: 
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g{xp) = g [r(x P ); p F , ie n ] = g [r; p F , is n ] 



(37) 



The commutator in the Eilenberger equations implies the existence of several 
invariants along the characteristic line r(x). For example, if the normalization 
condition Eq.(3) is fulfilled at a particular fixed point r(xo), it will be fulfilled 
everywhere along the line r(x). Likewise, the determinant detg(x) and the trace 
trg(x) remain constant for — oo < x < oo. 

Next we parametrize the 2x2 matrix Y(x) in the form 



Y = exp(a+K + ) exp(a 3 if 3 ) exp(a_if_ 



(38) 



in terms of three unknown functions a 3 (x) , a+(x) and a_(x) ( Eulcr like 'angles' 
in particle-hole space). The physical propagator, Eq.(36), assumes then the form 



9{x) 



-m ■ 



[1 - 2a_(x)a+(x)cxp(-a 3 (x))] • 2K 3 
+ a+(x) • [a_(x)a + (x) exp (— a 3 (x)) — 1] • 2K + 
+a_(x) • exp [— a 3 (x)] • 2K- 



(39) 



One finds from the differential equation for Y(x) a set of three coupled differential 
equations for a 3 (x) , a+(x) and a_(x) : 



a 3 — 2a + cxp(— a 3 ) a_ = — 

TlVF 

exp(-a 3 ) a_ = 

TlVF 

a + -a + a 3 +a, cxp(-a 3 ) a_ = 

TlVF 



(40) 
(41) 
(42) 



Here a (x) = ^a(x). It is readily seen that the three equations decouple, and 
that a_ and a 3 may be expressed in terms of a + only : 



a 3 (x) = - 



£ n X 



dsA^ (s)a + (s 
l /"^ 

a_(x) = — - — • / ds Z\'(s) exp [a 3 (s)] + a 



+4 0) 



(0) 



(43) 
(44) 



The differential equation that remains to be solved for a + (x) is a Riccati equa- 
tion: 

d 

hvFg-a + (x) + [2e n + A f (x) a+(x)] a+(x) - Z\(x) =0 (45) 

However, the accurate numerical calculation of the nested integral for a_(x) is 
time consuming (even on a fast computer). To overcome this difficulty we use a 
trick. 

Let 5a(x) and gs(x) be two different solutions of the Eilenberger equations. 
Then not only the linear combination ca gA (x) + cb ?b(x) is a solution , but the 
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products gs{x)- ()a(x) and ()a(x)- are solutions as well. For example, the 

linear combination gs (x) ■ gA (x) — gA {x) ■ gs (x) solves the Eilenberger equations 
and fulfills the necessary condition trg(x) = . 

Let us construct two particular zero trace solutions to the Eilenberger equa- 
tions: 



with 



g A {x) = Y A {x)-K_ ■ \Y A {x] 
g B (x)=Y B (x)-K + - \? B (x) 



(46) 
(47) 



Ya = cxp(a + K + ) exp(a 3 K 3 ) exp(a__ft'_) 
Y B = cxp(b-K-) cnp(b 3 K 3 ) exp(b+K + ) 



(48) 
(49) 



denoting two equivalent fundamental systems Y^(x)and Yb(x). The different 
order of factors in the defining expressions for Ya and Yb serves the purpose to 
avoid the difficult terms a_(x) and b+(x) in the expressions for gA and gB ■ The 
evaluation of nested integrals, see Eq.(43), is then not necessary. 

The set of equations fulfilled by b 3 (x) and b±(x) is only slightly different from 
the one for a 3 (x) and a±(x) : 



& 3 +26_cxp(fo 3 ) b+ = 

Tivf 

b- +b- & 3 +&_ exp(fo 3 ) b+ = 



(50) 
(51) 
(52) 



Here b (x) = -S^b{x). It is readily seen that the three equations decouple, and 
that 6+ (x) and 63 (x) may be expressed in terms of b- (x) only : 



n x + / dsA(s)b-(s) 
Jo 

b+{x) = ^—f dsA(s)e X p[-b 3 (s)] + b { ° 



(53) 
(54) 



The differential equation to be solved for 6-(x) is also a Riccati equation: 

hvF-^b-(x) - [2s n + A(x)b-(x)] b-{x) + A\x) = (55) 

We observe that any solution of this differential equation is related to the 
Riccati equation Eq.(55) via a reciprocity relation: 
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If a + (x) solves Eq.(45), then 

M * )= -^b (56) 

solves Eq.(55). 

We continue now our construction of the physical propagator. From the defin- 
ing equations, Eq.(46) and Eq.(47), we find the following explicit expressions: 

g A = cxp(-a 3 ) (K- - 2a+ K 3 + a\ K+^j (57) 

g B = exp(6 3 ) (K+ + 26_ K 3 + b 2 _ £_) (58) 

Note that the square of gA (x) and g B (x) vanishes identically, 

9A-9A = = g B ■ 9b , (59) 

because K± = . For x — > ±oo the propagators gA and gs 'explode', i.e. 



g A ,B~eM±^el + \A\ 2 ) (60) 

On the other hand, the commutator \gA(x),gB(x)] remains bounded in the limit 
x — > ±oo. The observation that a bounded solution to the Eilenberger equations 
may be constructed using the commutator of two unbounded solutions gA (x) and 
9b{x) is the well known 'explosion' trick [5]. 

The general (particle-hole symmetric) solution to the Eilenberger equations 
(2) may be given in the form 

g{x) = ca9a{x) + c B gB{x) + \§a(x), 9b{x)] 

Here, ca and cb represent initial values, 63(0) and 03(0) , to the functions b 3 (x) 
and a 3 (x). Of course, in an unbounded region exploding solutions must be for- 
bidden. Then the physical propagator g must be written, on either side of the 
turning point x — , as a superposition of a decaying solution and a bounded 
solution: 

7s(t\ - / Cb 9b(x) + \9a{x), g B {x)] if x > 

9[x) ~ \ ca9a{x) + [g A (x), g B (x)} if x < [bL > 

The square of g(x) is in this case independent on the constants ca and c B : 

g(x) -g(x) = [g A (x), g B (x)\ ■ \§a(x), gB{x)\ = -tt 2 ■ i (62) 

It is not difficult to show, that ca = = c B , provided the propagators 9a{x) 
and g B {x) are continuous at x — 0. 

In fact, let us assume the contrary: ca- c b 7^ . Continuity oig{x) at x = 
leads to 

cb9b(0 + ) + [?a(0+), ?b(0+)] - .9(0) - c A g A (0-) + [?a(0"), ?b(0")] (63) 
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Both solutions, gA(x) and g B {x) , are continuous at x — . Then it follows from 
Eq.(63) : c B g B (0) = ca<?a(0). This implies, in turn, a vanishing commutator, 
[<7a(0), 3b(0)] = 0, since gs(0) and 3a (0) become proportional. Also, the physical 
solution g(0) at x — must fulfill the normalization condition, i.e. g(0) • g(0) = 
-it 2 • f . But g B (0) -g B (0) = 6 = g A (0) ■ g A (Q) according to Eqs.(59). This is a 
contradiction! Hence Ca = = c B . 

The conclusion is that in an infinitely extended system the physical propaga- 
tor g(x) is completely determined by the commutator of the 'exploding' solutions: 



d{x) = [9A(x),g B (x)} = exp(6 3 - a 3 ) 



I-(a+fc_) 2 2ia + (l + a+b-) 
-2ib_ (1 + a+6_) -l + (a+fo_) 2 



Next we check the normalisation condition: 

g(x) -g(x) = [g A (x), g B (x)} ■ \§a(x), g B (x)} 
= \9i{x) ■ 33 (x) +g+ (x) -g-(x)] -T 
= [1 + a + (x)6_(x)] 4 • cxp[26 3 (x) - 2a 3 (x)] T 
= C- 1 



Indeed, C is a constant multiple of unity: 

d_ 

dx 



{[1 + a + (x)b_(x)} 4 ■ exp [26 3 (a;) - 2a 3 (x)}} = 



(64) 



(65) 



(66) 



From the normalisation condition, C = —n 2 , there follows for all x (up to a 
sign ± that is chosen to coincide with the bulk propagator) : 



exp [b 3 (x) - o 3 (a;)] = 



[1 + a + (x)b_(x)Y 
Then the Eilenberger propagator may be parametrised in the form: 



g(x) = 



1 + a(x) ■ b{x) 



1 — a(x) ■ b(x) 2i ■ a(x) 
-2i-b(x) -l + a(x)-b(x) 



(67) 



(68) 



Here and in the following we use notation such that &_ (x) = b{x) and a+(x) = 
a(x), since the other functions d-(x),a 3 (x) and b + (x),b 3 (x) are obsolete for the 
parametrisation of the Eilenberger propagator. It is remarkable that the solution 
to the Eilenberger equations (2) may be given a representation where it depends 
just on the solution of an initial value problem to a scalar differential equation 
of the Riccati type[15],[13]. 

To integrate the Riccati equations (45, 55) in a stable manner we need suitable 
initial values for the functions b(x) and a(x). For ie n situated in the upper half of 
the complex plane the function a(x) may be found in a stable manner integrating 
Eq.(45) as an initial value problem from x = — oo towards increasing a;- values, 
while the function b(x) may be found integrating Eq.(55) as an initial value 
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problem from x = +00 backwards towards decreasing z-values. The initial values 
for a(x) at x = —00 and b(x) at x = +00 are 

<,(_«,) = ^(-°°) (69) 

e n + ^sl + \A{-^)\ 2 

K+oo) = (70) 

l^(+oo)| 2 

provided ie n is in the upper half of the complex plane. 
The differential equations to be solved are: 

d 

Hv F —a(x) + [2e„ + A ] (x) ■ a(x)] ■ a(x) - A(x) = (71) 
Tiv F ^-b{x) - [2i n + A(x) ■ b(x)\ ■ b(x) + A\x) = (72) 

Sometimes knowledge of just one of the functions, say a(x), along a line r(x) 
(for —00 < x < 00) suffices to fix the other function, b(x), along the same line. 
An illustrative example is provided by a single cylindrically symmetric vortex 
line, orientated parallel to c, and centered at the origin of the ab -plane, say at 
R = 0. Due to energetic reasons, of course, only a single quantum of circulation, 
2^, is attached to the vortex. The corresponding pair potential becomes along 
the straight line r(x) = r a (x)a + ri,(x)h a function of x (and also of the impact 
parameter y) of the form: 

x + iy ie 



4(r( I ),p F )=F(7?T7,0). 



\/x 2 + y 2 



The prefactor F(^Jx 2 + y 2 , 9) is a suitable 'form factor' to shape the vortex core. 
We see from Eqs. (71,72) that in the presence of such a vortex line, b(x) is related 
to a(x) by symmetry: 

b(x) = -a(-x)e- 2lf> (73) 

Using Eq.(68) it follows that the corresponding Eilenberger propagator g for 
negative x is related to the propagator for positive x by the relation: 

g(-x) - -e^ • t 2 ■ g(x) ■ f 2 • e^ 3 (74) 

To determine the local density of states one needs the retarded and the 
advanced propagator of the quasiclassical theory. Actually, in equilibrium, only 
the retarded (or advanced) propagator is needed in the calculations, since both 
propagators are related to each other by complex conjugation. A convenient 
numerical method for the calculation of the retarded propagator g( re *) (r, 9, E) 
is to replace the discrete Matsubara frequency ie n according to the prescription 
ie n — > E + i0 + , and to solve the Riccati equations, Eqs. (71, 72), as functions of 
the energy E and of the impact parameters y and z. 
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If the denominator of the quasiclassical propagator, 1 + a(x) ■ b(x), becomes 
equal to zero at a point r(xo) for a characteristic energy E = Eb, it vanishes 
indeed for all x along the trajectory r(x): 



[1 + a(x) ■ b(x)] = exp 



-!- f ds (A(s)b(s)- A^s)a(s)) 
Tivf J x 



■ [1 + a(x Q ) ■ b(x )} 



(75) 

A simple proof of this relation uses differentiation with respect to x, and Eqs.(71, 
72). So, if the denominator 1 + a(x ) ■ b(x ) of the quasiclassical propagator 
Eq.(68), considered as a function of energy E, displays a simple zero at E = Eb, 
this zero, Eb, has a natural interpretation as a bound state energy, provided there 
exists a finite residue of the retarded propagator at E = Eb + i0 + . 

Since it is almost never possible to solve the equations of superconductivity 
exactly, one needs numerical methods. The alternative to a straightforward (but 
costly) numerical solution of the BdG-eigenvalue problem is the numerical solu- 
tion of the Eilenberger equations, provided the fundamental condition of quasi 
classical theory, Uf • £ 3> 1, is valid. For example, using the quasiclassical ap- 
proach of Eilenberger, it is comparatively easy to determine the deep lying bound 
states Eb of localised vortex core fermions (attached to a single vortex line) as 
a function of the impact parameters :Eb — Eb(y,z). The summation over (ex- 
act) eigenenergies of the bound states obtained from solving the BdG-eigenvalue 
problem (see Ref.[7]) for a single vortex line becomes equivalent to integrating 
the quasiclassical propagator with respect to the impact parameter. In certain 
materials, for example cuprates, the Cooper pairs display (perhaps) an uncon- 
ventional d x 2_ y 2 -symmetry. Results of a quasiclassical calculation of the bound 
state spectrum of quasiparticles around a single flux line in a superconductor 
with d x 2_ y 2 -pairing symmetry are published in Refs.[13],[14]. 

For superfluid 3 He — B, a prominent system with unconventional p-wave 
pairing symmetry, our method can be extended to the 4 x 4-Eilcnberger propa- 
gator for triplet pairing. A calculation of the spectrum of vortex core fermions 
around the o-vortex, the w-vortex and also the double core vortex is discussed 
in Rcfs. [16], [15]. 

We conclude that the quasiclassical propagator may be determined solving 
an initial value problem for a scalar differential equation of the Riccati type. For 
numerical calculations this method to solve the Eilenberger equations may be 
recomended for its intrinsic stability and speed. Also the Eilenberger approach 
is a suitable one for parallel computers. 



4 Connection to Linearized Bogoliubov-de Gennes 
Equations 

There exists an interesting connection between the solutions u(x) and v(x) to 
the linearized BdG-equations (often refered to as Andreev equations [3]), and 
the solutions to the Riccati equations. Along a characteristic line r(x) orientated 
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parallel to the Fermi velocity we have: 



-ihv F — 
Ox 



u{x) 




v{x) 





ie n (x) -A(x) 
A^{x) -is n (x) 



u(x) 
v(x) 



(76) 



We observe that a{x) , the solution to the Riccati Equation Eq.(45) , may be 
represented as the ratio: 

This connection between a(x) and the amplitudes u{x) and v{x) is readily 
demonstrated: 



-hvFTr-a(x) 
ox 



(78) 



v(x) [-ihv F -§^u{x)] -u(x) [-ihv F -^v(x)] 
v(x) \ie n {x)u(x) — A(x)v(x)} — u(x) [A^ (x)u(x) — ie n (x)v(x)] 

= [2i n + A\x) ■ a(x)] ■ a(x) - A(x) 

The conclusion is, that any solution to the Eilenberger equations (for ie n {x) 
in the upper half of the complex plane for x — > — oo ) may be reconstructed 
from a solution to the linearized BdG-diffcrential equations, provided the initial 
values for u(x) and v(x) are choosen in accordance with the known asymptotic 
behaviour of a(x) for x — > — oo (see Eq.(69)). 

Another interesting feature follows making a polar decomposition of the pair 
potential: 

A(x) = \A(x)\e i +W (79) 



and making the transformation 



a(x) 



in Eq.(45). It is readily shown that 
d 



v F —n{x) 
ox 



2ie n (x) + 2 \A{x)\ sin [t)(x) - 4>{x)\ = 



(80) 



(81) 



It is interesting that this differential equation, which is equivalent to Eq.(45), 
was already derived by J. Bardeen et al. [8] in their so called WKBJ-approach 
to the BdG-equations. 

So, the 'new' result is not the Riccati equation Eq.(45) itself (or any equiv- 
alent representation). New is the fact, that the full quasiclassical propagator 
g (r, 9, ie n ) (in equilibrium) is simply a rational function (see Eq.(68)) of the so- 
lutions to a scalar Riccati equation. In turn, solutions to these Riccati equations 
are related by Eq.(77) to the solutions of an initial value problem for the lin- 
earised Bogoiliubov-de Gennes equations. This implies, that the standard proce- 
dure of wave mechanics to calculate the observables of superconductivity, namely 
first solving the Bogolubov-de Gennes eigenvalue problem and afterwards doing 
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a summation over the eigenfunctions and eigenvalues to calculate the Green's 
function, is in deed obsolete and may be replaced by solving a scalar Riccati 
equation, or equivalently, by solving an initial value problem for the linearised 
BdG-equations, provided the quasiclassical condition k F £ 3> 1 is fulfilled. 



5 Exact Solutions 



In some cases we may calculate exact solutions to the Eilenberger equations for 
a given profile of the pair potential. First the term VfA(i) in the definition of 
£„(x) ( local Doppler shift) may be removed from the diagonal of Eqs.(76) mak- 
ing a gauge transformation, i.e. we may always assume i n (x) — > e n independent 
on x . Then, after a decomposition of the (transformed) pair potential into real 
and imaginary parts, 

A(x) = A x (x) + iA 2 (x) (82) 
the linearised BdG-equations (76) may be given the form: 



where 



d ^ ^ 

t 3 • ihv F — + A\{x) ■ ?i - A 2 (x) ■ ? 2 + is n ■ 1 



7, , ( u(x) 
\v(x) 



?sip{x) = 



(83) 
(84) 



Introducing an auxiliary spinor x(x) vm 

d 



T 3 ip(x) 



% ■ ^vft^; + Ai(x) ■ ?x - A 2 (x) ■ % - is n ■ 1 



X(x) (85) 



the following 2 nd order linear differential equation be derived: 

d 2 



2„,2 



n v 
—Hvf 



F dx 2 
(dAx{x) 



+ Ai(x) + Ai(x) + ei -1 



\ dx 



^ dA 2 {x) ^ 

t 2 H 5 — • n 

ox 



X(x) = 



(86) 



Provided an exact solution (in accordance with the initial conditions Eqs.(69), 
(70)) to Eq.(86) can be found, the quasiclassical propagator may be calculated 
exactly from Eqs. (77,68). This procedure is a useful principle for the construction 
of exact solutions to the Eilenberger equations. 

The major obstacle to decompose Eq.(86) just into two decoupled scalar 
differential equations is the spatial dephasing between imaginary and real parts 
of the gradient of the pair potential, 9A g^ and dA g^ , in Eq.(86). A simple 

case occurs, for example, if 9A q! c x ' > = . In this case a rotation around the fi-axis 
by an angle ^ leads to a diagonal matrix, i.e. the problem may be effectively 
decoupled into two scalar differential equations of 2 nd -order. The well known 
WKB method [9] in its standard guise may then be applied to construct an 
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(approximate) analytical solution in this case. If there is no dephasing between 
real and imaginary parts of the pair potential, some pair potentials A(x) with 
model character allow the construction of exact solutions, as we show below. 

On the other hand, if there exists dephasing between real and imaginary 
parts of the pair potential, the problem is more difficult. One way to proceed is 
stratification. This means one approximates the pair potential, A{x) = Ai(x) + 
i ■ A 2 (x), by a sequence of strata along x, such that A\(x) and A 2 (x) become 
continous and piecewise linear functions of x: 



A^x) 
A 2 (x) 



C\x + d\ 
c 2 x + d 2 



(87) 



The real constants d\,d 2 , c\,c 2 may change from one stratum to another stratum, 
like a (linear) spline function. Inside a fixed stratum the 2x2 matrix 



dMx) 
dx 



t 2 + 



dA 2 (x) 
dx 



Tl = ClT 2 + C 2 Tl 



(89) 



may be diagonalised using a suitable, x- independent unitary transformation ma- 
trix U of the form 



U 



1 

71 



T 3 + 



Cl 



V4 



T 2 + 



C 2 



ClT 2 



U Z = 1 

c 2 ?i = U ■ Jc\ +cl? 3 -U 



(90) 

(91) 
(92) 



Since the constants ci , c 2 may change from a given stratum to the neighbour- 
ing stratum, the unitary transformation matrix U changes accordingly. Inside a 
fixed stratum Eqs.(86) may be decomposed into two decoupled scalar differential 
equations of 2 nd -order for the components of the new spinor: 



U ■ x(x) = $(x) = 



(93) 



M x ). 

Then the transformed equations assume the familiar form of a 'shifted oscillator': 

d 2 



tfvp-^r + (ax + diY + (c 2 x + d 2 y + e) 



dx 2 



— tiVf 



<P(x) = 



(94) 



Exact solutions (inside a choosen stratum) may be presented to these differential 
equations using the linearly independent parabolic cylinder functions D „ (x) and 
D- u -i{ix). The latter special functions are solutions to the differential equation 
of the harmonic oscillator: 



dx 2 



1 

V+ 2 



D(x) 







(95) 



Transformation of the Eilenberger Equations 15 



The components of <&{x) are related to x{x) via the relation x{x) = U ■ <P(x) , 
i.e. the changes of the constants c\ and C2 from one stratum to the next stratum 
determine also the admixture of <&i{x) and #2 0*0 , and henceforth the spinor 

X(x): 



V2 



X2(X) = T2 



V c i + 4 

C2 + ici 

V c i + 4 



$i{x) 



(96) 
(97) 



Finally, from Eq.(85) the ratio a(x) — i ■ may be calculated such that a (a;) 
remains a smooth function of x . Once b(x) has been found from a similar con- 
sideration (or from a(x) by a symmetry argument) the quasiclassical propagator 
follows from Eq.(68). 

Closing this section we give three examples of pair potential profiles A(x) 
that allow an exact solution to the Eilenberger equations. Allthough these pair 
potential are not selfconsistent they are helpful for a qualitative physical under- 
standing: 



A(x) = A c 
A(x) = A c 



x + iy 



iy 



\/x 2 + y 2 



A(x) = A x ■ tanh(-) 



(98) 
(99) 
(100) 



The model Eq.(98) represents the inner core of a vortex, and an exact solu- 
tion to the Eilenberger equations may be constructed in terms of the parabolic 
cylinder functions along the lines explained above [15]. 

The model Eq. (99) represents the outer core of a vortex ( it describes a 
pure phase vortex), and may actually be solved exactly for the special case 
ie n -> E = IAJ. 



a(x) — 



(1 - 2iW) (y + ^x 2 + y 2 ) + (1 + 2iW) ■ ix 
(1 + 2iW) (y + ^x 2 + y 2 ) - (1 - 2iW) ■ ix 



I A, 



(101) 

(102) 
(103) 



At the gap edge, E = |A»| > the corresponding expression for the quasiclassical 
propagator, Eq. (68), displays algebraic decay for x — > ±00. Certainly, it would 
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be nice to know the exact solution also for other energies E , but this problem 
is still unsolved. Nevertheless, a qualitatively correct (but not exact) solution 
to Eq.(86) for the case of a vortex line may be found [15] using a method of 
Stueckelberg[10]. 

The third model, Eq.(lOO), is simpler than the previous model, Eq.(99), be- 
cause the pair potential A(x) is a real function (no dephasing of components of 
spinor x). The differential equation for x{x), Eq.(86), may be solved exactly in 
terms of the hypergeomctric function [19]). A particularly simple looking analytic 
solution results if the size parameter £ of the domain wall described by Eq.(lOO) 
assumes the special value £ = ^p 1 : 

e n - y/el + Al, + A x tanhf 

a(x) = 7 2 I A2 a I ( 104 ) 

£n - V4 + A 2 OQ -A oc tanh | 



A 2 



2e n cosh 2 



T3 



A 



3o tanh | • t 2 



' 2e„ 



(105) 



We see that the explicit analytical expression for the quasiclassical propagator, 
g(x) , nicely reveals the existence of a mid gap state (a quasiparticle bound state 
with an excitation energy E around zero), assuming we make the analytical 
continuation ie n — ► E+i0 + . The existence of such a mid gap state is in agreement 
with our previous remarks associated with Eq.(75). 



6 Outlook 

In this article we have shown how to parametrise the 2 x 2-Eilcnbcrgcr equations 
in terms of the solutions to a scalar Riccati equation, or equivalcntly, in terms 
of the solutions to an initial value problem for the linearised BdG-equations. 
The latter method for the reconstruction of the Eilenberger propagator may be 
extended [15] to the full 4 x 4-matrix equations of the quasiclassical theory of 
superconductivity and superfluidity (including paramagnetic effects). For a gen- 
eralisation of the method to non equilibrium see [17], [20]). For an authoritative 
review of the quasiclassical theory see Ref . [4] . 
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